home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / CUGUK / APPLICAT / C050.ZIP / SFSSRC.ZIP / SFS.SRC / AS / AS_ORBIT.C < prev    next >
C/C++ Source or Header  |  1991-11-04  |  10KB  |  365 lines
  1. /****************************************************************
  2.  
  3.     as_orbit.c      Orbital calculation routines
  4.             for astronomy (as) subsystem
  5.  
  6.             Copyright (c) 1991, Ted A. Campbell
  7.  
  8.             Bywater Software
  9.             P. O. Box 4023 
  10.             Duke Station 
  11.             Durham, NC  27706
  12.  
  13.             email: tcamp@hercules.acpub.duke.edu
  14.  
  15.     Copyright and Permissions Information:
  16.  
  17.     All U.S. and international copyrights are claimed by the
  18.     author. The author grants permission to use this code
  19.     and software based on it under the following conditions:
  20.     (a) in general, the code and software based upon it may be 
  21.     used by individuals and by non-profit organizations; (b) it
  22.     may also be utilized by governmental agencies in any country,
  23.     with the exception of military agencies; (c) the code and/or
  24.     software based upon it may not be sold for a profit without
  25.     an explicit and specific permission from the author, except
  26.     that a minimal fee may be charged for media on which it is
  27.     copied, and for copying and handling; (d) the code must be 
  28.     distributed in the form in which it has been released by the
  29.     author; and (e) the code and software based upon it may not 
  30.     be used for illegal activities. 
  31.  
  32. ****************************************************************/
  33.  
  34. #include "math.h"
  35. #include "time.h"
  36. #include "bw.h"
  37. #include "as.h"
  38.  
  39. #define KEP_ACCURACY    0.00001
  40.  
  41. /****************************************************************
  42.  
  43.    or_init()      Set up an orbit
  44.  
  45. ****************************************************************/
  46.  
  47. or_init( orbit, ap, aa )
  48.    struct as_orbit *orbit;
  49.    double aa, ap;
  50.    {
  51.  
  52. #ifdef  DEBUG
  53.    if ( orbit->focus->radius <= 0.0 )
  54.       {
  55.       sprintf( bw_ebuf, "[pr:] or_init() received orbit->focus->radius = %lf",
  56.      orbit->focus->radius );
  57.       bw_error( bw_ebuf );
  58.       return;
  59.       }
  60.    if ( orbit->focus->sid_day < 0.0 )
  61.       {
  62.       sprintf( bw_ebuf, "[pr:] or_init() received orbit->focus->sid_day = %lf",
  63.      orbit->focus->sid_day );
  64.       bw_error( bw_ebuf );
  65.       return;
  66.       }
  67.    if ( orbit->focus->mass < 0.0 )
  68.       {
  69.       sprintf( bw_ebuf, "[pr:] or_init() received orbit->focus->mass = %lf",
  70.      orbit->focus->mass );
  71.       bw_error( bw_ebuf );
  72.       return;
  73.       }
  74.    if ( ( orbit->inclination > 180.0 ) || ( orbit->inclination < -180.0 ))
  75.       {
  76.       sprintf( bw_ebuf, "[pr:] or_init() received orbit->inclination = %lf",
  77.      orbit->inclination );
  78.       bw_error( bw_ebuf );
  79.       return;
  80.       }
  81.    if ( ( orbit->lon_an > 180.0 ) || ( orbit->lon_an < -180.0 ))
  82.       {
  83.       sprintf( bw_ebuf, "[pr:] or_init() received orbit->lon_an = %lf",
  84.      orbit->lon_an );
  85.       bw_error( bw_ebuf );
  86.       return;
  87.       }
  88.    if ( ( orbit->arg_per > 180.0 ) || ( orbit->arg_per < -180.0 ))
  89.       {
  90.       sprintf( bw_ebuf, "[pr:] or_init() received orbit->arg_per = %lf",
  91.      orbit->arg_per );
  92.       bw_error( bw_ebuf );
  93.       return;
  94.       }
  95. #endif
  96.  
  97.    /***    Calculate semimajor axis        */
  98.  
  99.    orbit->semimajor = ( aa + ap + ( 2 * orbit->focus->radius )) / 2;
  100.  
  101.    /***    Calculate orbital center - focus distance       */
  102.  
  103.    orbit->dist = orbit->semimajor - ( orbit->focus->radius + ap );
  104.  
  105.    /***    Calculate semiminor axis        */
  106.  
  107.    orbit->semiminor = sqrt( (double) ( orbit->semimajor * orbit->semimajor ) - (double) ( orbit->dist * orbit->dist ));
  108.  
  109.    /***    Calculate orbital eccentricity  */
  110.  
  111.    orbit->eccentricity = sqrt( 1.0 - (( orbit->semiminor / (double) orbit->semimajor ) * ( orbit->semiminor / (double) orbit->semimajor )) );
  112.  
  113.    /***    Calculate orbital period        */
  114.  
  115.    orbit->period = sqrt ( ( ( 4.0 * ( PI * (double) PI ) ) / ( UGC * orbit->focus->mass ) )
  116.       * ( orbit->semimajor * orbit->semimajor * orbit->semimajor ) );
  117.  
  118.    /***   Calculate the increment factor of longitude of 
  119.       the ascending node.  This factor must be multiplied
  120.       by a time factor (orbital period or time into orbit, 
  121.       in seconds).  See Davidoff, pp. 8-9 and 8-10, and
  122.       formula 8.14.    */
  123.  
  124.    if ( orbit->focus->sid_day == 0.0 )
  125.       {
  126.       orbit->lif = 0;
  127.       }
  128.    else
  129.       {
  130.       orbit->lif = (2 * PI) / (double) orbit->focus->sid_day;
  131.       }
  132.  
  133.    }
  134.  
  135. /****************************************************************
  136.  
  137.    or_kep()      Solve Kepler's equation
  138.  
  139.    Globals utilized:  
  140.  
  141.    orbit->eccentricity Eccentricity of the orbital ellipse
  142.  
  143.    orbit->semimajor    Semimajor axis of the orbital ellipse (km)
  144.  
  145.    orbit->period       Orbital period (seconds)
  146.  
  147.    Input:
  148.  
  149.    t                   Time from periapsis in seconds
  150.                ( 0 < t < orbit->period )
  151.  
  152.    Output:
  153.  
  154.    theta               Angle between periapsis, geocenter, and
  155.                current position (theta).
  156.  
  157.    r                   Distance from satellite to focal center,
  158.                in kilometers
  159.  
  160. ****************************************************************/
  161.  
  162. or_kep( orbit, t, theta, r )
  163.    struct as_orbit *orbit;
  164.    long t;
  165.    double *theta;
  166.    long *r;
  167.    {
  168.    double z, z3, E;
  169.    z  = 2.0 * PI * ( (double) t / (double) orbit->period );
  170.    E  = z;
  171.  
  172.    do
  173.       {
  174.       z3 = ( E - ( orbit->eccentricity * sin( E )) - z ) / ( 1 - ( orbit->eccentricity * cos( E )));
  175.       E  = E - z3;
  176.       }
  177.    while ( fabs( z3 ) > KEP_ACCURACY );
  178.  
  179.    *theta = PI;
  180.    if ( E != PI )
  181.       {
  182.       *theta = 2.0 * atan( sqrt( ( 1 - orbit->eccentricity ) / ( 1 + orbit->eccentricity )) * tan( E / 2.0 ));
  183.       if ( E > PI )
  184.          {
  185.          *theta = ( (double) 2.0 * PI ) + *theta;
  186.          }
  187.       }
  188.  
  189.    *r = orbit->semimajor * ( 1 - orbit->eccentricity * orbit->eccentricity ) / ( 1 + orbit->eccentricity * cos( *theta ));
  190.    return 1;
  191.    }
  192.  
  193.  
  194. /****************************************************************
  195.  
  196.    or_ssp()   Calculate Subsatellite Point
  197.  
  198.    This function utilizes available data on the satellite 
  199.    to return the subsatellite point, that is, the point 
  200.    on the surface of the orbital focus directly under 
  201.    the satellite at a given moment.  
  202.  
  203. ****************************************************************/
  204.  
  205. or_ssp( orbit, t, latitude, longitude, r, n )
  206.    struct as_orbit *orbit;
  207.    long t, *r, *n;
  208.    double *longitude, *latitude;
  209.    {
  210.    static long _r, tp;
  211.    static double theta;
  212.    double n1, n2, n4, E, lan;
  213.    long t1;
  214.  
  215. #ifdef  DEBUG
  216.    if ( orbit->focus->radius < OR_RAD_MIN )
  217.       {
  218.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->focus->radius = %lf",
  219.      orbit->focus->radius );
  220.       bw_error( bw_ebuf );
  221.       return;
  222.       }
  223.    if ( orbit->focus->sid_day < OR_SID_MIN )
  224.       {
  225.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->focus->sid_day = %lf",
  226.      orbit->focus->sid_day );
  227.       bw_error( bw_ebuf );
  228.       return;
  229.       }
  230.    if ( orbit->focus->mass < OR_MAS_MIN )
  231.       {
  232.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->focus->mass = %le",
  233.      orbit->focus->mass );
  234.       bw_error( bw_ebuf );
  235.       return;
  236.       }
  237.    if ( ( orbit->inclination > OR_INC_MAX ) || ( orbit->inclination < OR_INC_MIN ))
  238.       {
  239.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->inclination = %lf",
  240.      orbit->inclination );
  241.       bw_error( bw_ebuf );
  242.       return;
  243.       }
  244.    if ( ( orbit->lon_an > OR_LAN_MAX ) || ( orbit->lon_an < OR_LAN_MIN ))
  245.       {
  246.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->lon_an = %lf",
  247.      orbit->lon_an );
  248.       bw_error( bw_ebuf );
  249.       return;
  250.       }
  251.    if ( ( orbit->arg_per > OR_ARP_MAX ) || ( orbit->arg_per < OR_ARP_MIN ))
  252.       {
  253.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->arg_per = %lf",
  254.      orbit->arg_per );
  255.       bw_error( bw_ebuf );
  256.       return;
  257.       }
  258.    if ( orbit->period < OR_PER_MIN )
  259.       {
  260.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->period = %lf",
  261.      orbit->period );
  262.       bw_error( bw_ebuf );
  263.       return;
  264.       }
  265.    if ( ( orbit->eccentricity < OR_ECC_MIN ) || ( orbit->eccentricity > OR_ECC_MAX ))
  266.       {
  267.       sprintf( bw_ebuf, "[pr:] or_ssp() received orbit->eccentricity = %lf",
  268.      orbit->eccentricity );
  269.       bw_error( bw_ebuf );
  270.       return;
  271.       }
  272. #endif
  273.  
  274.    *n = t / (long) orbit->period;           /* return number of current orbit */
  275.    t1 = t - ( orbit->period * ( *n ) );     /* get seconds into this orbit */
  276.    if ( t1 == 0 )
  277.       {
  278.       t1 = 1;
  279.       }
  280.  
  281.    or_kep( orbit, t1, &theta, &_r );
  282.    *r = _r;
  283.  
  284. #ifdef  OLD_DEBUG
  285.     sprintf( bw_ebuf, "or_ssp(): r     = %ld (%ld) ", _r, *r );
  286.     bw_debug( bw_ebuf );
  287.     sprintf( bw_ebuf, "or_ssp(): theta = %lf", theta );
  288.     bw_debug( bw_ebuf );
  289. #endif
  290.  
  291.    /***   Calculate the longitude of the ascending node at 
  292.       the beginning of this orbit.  See Davidoff, p. 8-9
  293.       and 8-10, and equations 8.13a, 8.13b, and 8.14.   */
  294.  
  295.    lan = ( orbit->lon_an * DEG_RAD ) - ( orbit->lif * orbit->period * *n );
  296.  
  297.    /***   Calculate the latitude of the SSP.  See Davidoff, 
  298.       pp. 8-13 - 8-15, esp. equation 8.20.      */
  299.  
  300.    *latitude = asin( sin( orbit->inclination * DEG_RAD )
  301.       * sin( theta + ( orbit->arg_per * DEG_RAD )));
  302.  
  303.    if ( ( orbit->inclination >= 0 ) && ( orbit->inclination < 90 ) )
  304.       {
  305.       n2 = 1;
  306.       }
  307.    else
  308.       {
  309.       n2 = 0;
  310.       }
  311.    if ( ( ( *latitude ) * RAD_DEG ) < 0.0 )
  312.       {
  313.       n4 = 1;
  314.       }
  315.    else
  316.       {
  317.       n4 = 0;
  318.       }
  319.  
  320.    if ( ( orbit->arg_per > 180 ) && ( orbit->arg_per <= 540 ) )
  321.       {
  322.       n1 = 1;
  323.       }
  324.    else
  325.       {
  326.       n1 = 0;
  327.       }
  328.  
  329.    E = 2 * atan( (double) pow( ( 1 - orbit->eccentricity ) / ( 1 + orbit->eccentricity ), (double) 0.5 )
  330.         * tan( ( orbit->arg_per * DEG_RAD ) / 2.0 )) + ( 2.0 * PI * n1);
  331.    tp = ( orbit->period / ( 2.0 * PI ) ) * ( E - sin( E ) );
  332.  
  333.    *longitude = as_lfix( lan
  334.            - pow( (double) -1.0, n2 + n4 )
  335.        * acos( cos( theta + ( orbit->arg_per * DEG_RAD ))
  336.        / cos( *latitude ) )
  337.        - orbit->lif * (double) t1
  338.        - orbit->lif * (double) tp );
  339.  
  340.    /* convert to degrees */
  341.  
  342.    *longitude = *longitude * RAD_DEG;
  343.    *latitude  = *latitude  * RAD_DEG;
  344.  
  345.    }
  346.  
  347. double
  348. as_lfix( l )
  349.    double l;
  350.    {
  351.    double l1;
  352.    l1 = l;
  353.    while ( l1 < ( 0 - PI ) )
  354.       {
  355.       l1 += ( 2.0 * PI );
  356.       }
  357.    while ( l1 > PI )
  358.       {
  359.       l1 -= ( 2.0 * PI );
  360.       }
  361.    return l1;
  362.    }
  363.  
  364.  
  365.